# R code to replicate the map generation
library(foreign)
library(cshapes)
library(maptools)
library(classInt)

# This program requires the Fearon & Laitin (2003) and Kristian Gleditsch's
# updated GDP data to be present in the working directory. These files can be
# obtained from
# http://www.stanford.edu/~jfearon/data/apsr03repdata.zip and
# http://ibs.colorado.edu/~ksg/trade/exptradegdpv4.1.zip (pick file "pwt_cow.asc" from the archive)

# modify this line, enter your working directory
setwd("/Users/nilsw/Projects/CShapes/R")

# load Fearon&Laitin (2003) dataset
fldata <- read.dta("repdata.dta", convert.factors=F)
fldata$onset[fldata$onset>1] <- 1 # fix coding error
fldata$ccode[fldata$ccode==260 & fldata$year>=1990] <- 255 # fix wrong code for Germany after 1990
fldata <- fldata[fldata$year>=1950,]
fldata$uniqueid <- fldata$ccode*10000 + fldata$year

# gdpenl has many missings. Replace with KSG trade and gdp data.
ksgdata <- read.table("pwt_cow.asc", header=T) 
ksgdata$uniqueid <- ksgdata$statenum*10000 + ksgdata$year
ksgdata <- subset(ksgdata, select=c(uniqueid, gdppc))


fldata <- merge(fldata, ksgdata, by="uniqueid")

fldata <- subset(fldata, 
  !is.na(onset) 
  & !is.na(warl) 
  & !is.na(gdppc) 
  & !is.na(lpopl1) 
  & !is.na(lmtnest) 
  & !is.na(ncontig) 
  & !is.na(Oil)
  & !is.na(nwstate)
  & !is.na(instab)
  & !is.na(polity2l)
  & !is.na(ethfrac)
  & !is.na(relfrac))
  
# estimate model 1
m1 <- glm(onset ~ warl+log10(gdppc)+lpopl1+lmtnest+ncontig+Oil+nwstate+instab+polity2l+ethfrac+relfrac, family=binomial, data=fldata)
fldata$yhat <- predict(m1, type = "response")
fldata.1985 <- fldata[fldata$year==1985,]
fldata.1995 <- fldata[fldata$year==1995,]

# prepare cshapes data frames for 1985 and 1995
cshp <- cshp()
startdate <- as.Date(paste(cshp$COWSYEAR, cshp$COWSMONTH, cshp$COWSDAY, sep="-"))
enddate <- as.Date(paste(cshp$COWEYEAR, cshp$COWEMONTH, cshp$COWEDAY, sep="-"))
date.1 <- as.Date("1985-6-30")
cshp.1985 <- cshp[startdate <= date.1 & enddate> date.1,]
date.2 <- as.Date("1995-6-30")
cshp.1995 <- cshp[startdate <= date.2 & enddate> date.2,]

# find cases that are present in both datasets
cshp.1985 <- cshp.1985[cshp.1985$COWCODE %in% fldata.1985$ccode,]
fldata.1985 <- fldata.1985[fldata.1985$ccode %in% cshp.1985$COWCODE,]
cshp.1995 <- cshp.1995[cshp.1995$COWCODE %in% fldata.1995$ccode,]
fldata.1995 <- fldata.1995[fldata.1995$ccode %in% cshp.1995$COWCODE,]

# join data 
o <- match(cshp.1985$COWCODE, fldata.1985$ccode)
fldata.1985 <- fldata.1985[o,]
cshp.1985 <- spCbind(cshp.1985, fldata.1985$yhat)
o <- match(cshp.1995$COWCODE, fldata.1995$ccode)
fldata.1995 <- fldata.1995[o,]
cshp.1995 <- spCbind(cshp.1995, fldata.1995$yhat)

# plot
pal <- grey.colors(4, 1, 0, 2.2)
q5 <- classIntervals(cshp.1985$fldata.1985.yhat, n=5, style="equal")
pdf("flpredictions1985.pdf")
plot(cshp.1985, col=findColours(q5, pal))
dev.off()
q5 <- classIntervals(cshp.1995$fldata.1995.yhat, n=5, style="equal")
pdf("flpredictions1995.pdf")
plot(cshp.1995, col=findColours(q5, pal))
dev.off()




